In [1]:
def odelay(func, y0, xspan, events, TOLERANCE = 1e-6, fsolve_args=None, **kwargs):
'''Solve an ODE with events.
func is callable, with signature func(Y, x)
y0 are the initial conditions xspan is what you want to integrate
over
events is a list of callable functions with signature event(Y, x).
These functions return zero when an event has happened.
TOLERANCE is what is used to identify when an event has occurred.
[value, isterminal, direction] = event(Y, x)
value is the value of the event function. When value = 0, an event
is triggered
isterminal = True if the integration is to terminate at a zero of
this event function, otherwise, False.
direction = 0 if all zeros are to be located (the default), +1
if only zeros where the event function is increasing, and -1 if
only zeros where the event function is decreasing.
fsolve_args is a dictionary of options for fsolve
kwargs are any additional options you want to send to odeint.
Returns [x, y, te, ye, ie]
x is the independent variable array
y is the solution
te is an array of independent variable values where events occurred
ye is an array of the solution at the points where events occurred
ie is an array of indices indicating which event function occurred.
'''
if 'full_output' in kwargs:
raise Exception('full_output not supported as an option')
if fsolve_args is None:
fsolve_args = {}
x0 = xspan[0] # initial point
X = [x0]
sol = [y0]
TE, YE, IE = [], [], [] # to store where events occur
# initial value of events
e = np.zeros((len(events), len(xspan)))
for i,event in enumerate(events):
e[i,0], isterminal, direction = event(y0, x0)
# now we step through the integration
for i, x1 in enumerate(xspan[0:-1]):
x2 = xspan[i + 1]
f1 = sol[i]
f2 = odeint(func, f1, [x1, x2], **kwargs)
X += [x2]
sol += [f2[-1,:]]
# check event functions. At each step we compute the event
# functions, and check if they have changed sign since the
# last step. If they changed sign, it implies a zero was
# crossed.
for j, event in enumerate(events):
e[j, i + 1], isterminal, direction = event(sol[i + 1], X[i + 1])
if ((e[j, i + 1] * e[j, i] < 0) # sign change in
# event means zero
# crossing
or np.abs(e[j, i + 1]) < TOLERANCE # this point is
# practically 0
or np.abs(e[j, i]) < TOLERANCE):
xLt = X[-1] # Last point
fLt = sol[-1]
eLt = e[j, i+1]
# we need to find a value of x that makes the event zero
def objective(x):
# evaluate ode from xLT to x
txspan = [xLt, x]
tempsol = odeint(func, fLt, txspan, **kwargs)
sol = tempsol[-1, :]
val, isterminal, direction = event(sol, x)
return val
from scipy.optimize import fsolve
xZ, = fsolve(objective, xLt, **fsolve_args) # this should be the
# value of x that makes
# the event zero
# now evaluate solution at this point, so we can
# record the function values here.
txspan = [xLt, xZ]
tempsol = odeint(func, fLt, txspan, **kwargs)
fZ = tempsol[-1,:]
vZ, isterminal, direction = event(fZ, xZ)
COLLECTEVENT = False
if direction == 0:
COLLECTEVENT = True
elif (e[j, i + 1] > e[j, i] ) and direction == 1:
COLLECTEVENT = True
elif (e[j, i + 1] < e[j, i] ) and direction == -1:
COLLECTEVENT = True
if COLLECTEVENT:
TE.append(xZ)
YE.append(fZ)
IE.append(j)
if isterminal:
X[-1] = xZ
sol[-1] = fZ
return (np.array(X),
np.array(sol),
np.array(TE),
np.array(YE),
np.array(IE))
# at the end, return what we have
return (np.array(X),
np.array(sol),
np.array(TE),
np.array(YE),
np.array(IE))
In [5]:
from pycse import *
import numpy as np
def myode(f, x):
return 3*x**2 + 12*x -4
def event1(f, x):
'an event is when f = 0 and event is decreasing'
isterminal = True
direction = -1
return f, isterminal, direction
def event2(f, x):
'an event is when f = 0 and increasing'
isterminal = False
direction = 1
return f, isterminal, direction
f0 = -120
xspan = np.linspace(-8, 4)
X, F, TE, YE, IE = odelay(myode, f0, xspan, events=[event1, event2])
import matplotlib.pyplot as plt
plt.plot(X, F, '.-')
# plot the event locations.use a different color for each event
colors = 'rg'
for x,y,i in zip(TE, YE, IE):
plt.plot([x], [y], 'o', color=colors[i])
plt.savefig('images/event-ode-2.png')
plt.show()
print TE, YE, IE
In [6]:
import numpy as np
from scipy.integrate import odeint
def dCadt(Ca, t):
"the ode function"
k = 0.23
return -k * Ca**2
Ca0 = 2.3
# create lists to store time span and solution
tspan = [0, ]
sol = [Ca0,]
i = 0
while i < 100: # take max of 100 steps
t1 = tspan[i]
Ca = sol[i]
# pick the next time using a Newton-Raphson method
# we want f(t, Ca) = (Ca(t) - 1)**2 = 0
# df/dt = df/dCa dCa/dt
# = 2*(Ca - 1) * dCadt
t2 = t1 - (Ca - 1.0)**2 / (2 * (Ca - 1) *dCadt(Ca, t1))
f = odeint(dCadt, Ca, [t1, t2])
if np.abs(Ca - 1.0) <= 1e-4:
print 'Solution reached at i = {0}'.format(i)
break
tspan += [t2]
sol.append(f[-1][0])
i += 1
print 'At t={0:1.2f} Ca = {1:1.3f}'.format(tspan[-1], sol[-1])
import matplotlib.pyplot as plt
plt.plot(tspan, sol, 'bo')
plt.show()
In [ ]: